DrainageDensity Function

public function DrainageDensity(channel, fdir, x, y, flowAcc) result(dd)

Drainage Density (1/m) of a basin is the total line length of all streams divided by basin area. If flow accumulation is not given, it is computed from flow direction. If coordinate of closing section is not given drainage density is computed for the entire grid.

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: channel

mask of channel cells. NoData for non channel cells
type(grid_integer), intent(in) :: fdir

flow direction
real(kind=float), intent(in), optional :: x

coordinate of basin closing section

real(kind=float), intent(in), optional :: y

coordinate of basin closing section

type(grid_real), intent(in), optional :: flowAcc

flow accumulation

Return Value real(kind=float)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: cLength

cell length
logical, public :: check
type(grid_real), public :: facc
integer(kind=short), public :: i
integer(kind=short), public :: is

coordinate of downstream cell in local reference system
integer(kind=short), public :: j
integer(kind=short), public :: js

coordinate of downstream cell in local reference system
real(kind=float), public :: length

total river length
type(grid_integer), public :: mask

Source Code

FUNCTION DrainageDensity &
!
(channel, fdir, x, y, flowAcc) &
!
RESULT (dd)

IMPLICIT NONE

!Arguments with intent (in)
TYPE(grid_integer),INTENT(IN) :: channel !!mask of channel cells. NoData for non channel cells
TYPE(grid_integer),INTENT(IN) :: fdir !!flow direction
REAL (KIND = float), INTENT(IN),OPTIONAL :: x, y !!coordinate of basin closing section 
TYPE (grid_real),INTENT(IN),OPTIONAL :: flowAcc !!flow accumulation

!local declarations
REAL (KIND = float) :: dd
TYPE (grid_real) :: facc
TYPE (grid_integer) :: mask
INTEGER (KIND = short) :: i,j
INTEGER (KIND = short) :: is, js !!coordinate of downstream cell in local reference system
REAL (KIND = float) :: length !!total river length
REAL (KIND = float) :: cLength !!cell length
LOGICAL :: check

!------------------------------end of declaration -----------------------------
!set flowaccumulation
IF (PRESENT (flowAcc)) THEN
   CALL NewGrid (facc, fdir)
   facc = flowAcc

ELSE
  CALL FlowAccumulation (fdir, facc)
END IF

IF (PRESENT(x) .AND. PRESENT(y)) THEN !delineate watershed mask
   CALL BasinDelineate (fdir,x,y, mask)
ELSE  !drainage density is computed on entire grid
   CALL NewGrid (mask, fdir)
END IF

!compute total channel length [m]
length = 0.
DO i = 1, mask % idim
  DO j = 1, mask % jdim
    IF (mask % mat (i,j) /= mask % nodata) THEN
       IF (channel % mat (i,j) /= channel % nodata) THEN
          CALL DownstreamCell (i, j, fdir % mat(i,j), is, js, dx = cLength, grid = fdir)
          length = length + cLength
       END IF
    END IF
  END DO
END DO

!compute drainage density [1/m]
CALL GetIJ (x, y, fdir, i, j, check)

IF (.NOT. check) THEN
  CALL Catch ('error', 'Morphology', 'Drainage density computation: ', &
             argument = 'closing section outside grid dimension')
END IF

IF (channel % mat(i,j) == channel % nodata) THEN
  CALL Catch ('error', 'Morphology', 'Drainage density computation: ', &
             argument = 'closing section on hillslope')
END IF

dd = length / facc % mat(i,j) 

!deallocate
CALL GridDestroy (facc)
CALL GridDestroy (mask)

RETURN

END FUNCTION DrainageDensity